home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / kw_test.pro < prev    next >
Text File  |  1997-07-08  |  6KB  |  167 lines
  1. ;$Id: kw_test.pro,v 1.4 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       KW_TEST
  8. ;
  9. ; PURPOSE:
  10. ;       This function tests the hypothesis that three or more sample 
  11. ;       popultions have the same mean of distribution against the
  12. ;       hypothesis that they differ. The popultions may be of equal
  13. ;       or unequal lengths. The result is a two-element vector containing 
  14. ;       the test statistic H and the one-tailed probability of obtaining 
  15. ;       a value of H or greater from a chi-square distribution. This type 
  16. ;       of test is often refered to as the Kruskal-Wallis H-Test.
  17. ;
  18. ; CATEGORY:
  19. ;       Statistics.
  20. ;
  21. ; CALLING SEQUENCE:
  22. ;       Result = KW_test(X)
  23. ;
  24. ; INPUTS:
  25. ;       X:    An array of m-columns (m >= 3) and n-rows of type integer, 
  26. ;             float or double. The columns of this two dimensional array 
  27. ;             correspond to the sample popultions. If the sample popultions
  28. ;             are of unequal length, all vectors must be appended up to a
  29. ;             common length of n using a user-specified missing value. This
  30. ;             method requires the use of the MISSING keyword.
  31. ;
  32. ; KEYWORD PARAMETERS:
  33. ;      DF:    Use this keyword to specify a named variable which returns
  34. ;             the degrees of freedom used to compute the probability of
  35. ;             obtaining a value of H or greater from the corresponding 
  36. ;             chi-square distribution
  37. ;
  38. ; MISSING:    Use this keyword to specify a non-zero numeric value which
  39. ;             is used to appended popultions of unequal length up to a 
  40. ;             common length of n.
  41. ;
  42. ; EXAMPLE:
  43. ;       Test the hypothesis that three sample popultions have the same mean 
  44. ;       of distribution against the hypothesis that they differ at the 0.05 
  45. ;       significance level.
  46. ;         sp0 = [24.0, 16.7, 22.8, 19.8, 18.9]
  47. ;         sp1 = [23.2, 19.8, 18.1, 17.6, 20.2, 17.8]
  48. ;         sp2 = [18.2, 19.1, 17.3, 17.3, 19.7, 18.9, 18.8, 19.3]
  49. ;       Since the sample popultions are of unequal lengths, a missing value
  50. ;       must be appended to sp0 and sp1. In this example the missing value
  51. ;       is -1.0 and the 3-column, 8-row input array is defined as:
  52. ;         x = [[24.0, 23.2, 18.2], $
  53. ;              [16.7, 19.8, 19.1], $
  54. ;              [22.8, 18.1, 17.3], $
  55. ;              [19.8, 17.6, 17.3], $
  56. ;              [18.9, 20.2, 19.7], $
  57. ;              [-1.0, 17.8, 18.9], $
  58. ;              [-1.0, -1.0, 18.8], $
  59. ;              [-1.0, -1.0, 19.3]]
  60. ;         result = kw_test(x, missing = -1)
  61. ;       The result should be the 2-element vector:
  62. ;         [1.65862, 0.436351]
  63. ;       The computed probability (0.436351) is greater than the 0.05
  64. ;       significance level and therefore we do not reject the hypothesis
  65. ;       that the three sample popultions s0, s1 and s2 have the same mean 
  66. ;       of distribution.
  67. ;
  68. ; PROCEDURE:
  69. ;       KW_TEST computes the nonparametric Kruskal-Wallis H-Test for three or
  70. ;       more populations of equal or unequal size. This test is an extension
  71. ;       of the Rank Sum Test implemented in the RS_TEST function. When each 
  72. ;       sample population contains at least five observations, the H test
  73. ;       statistic is approximated very well by a chi-square distribution with
  74. ;       DF degrees of freedom. The hypothesis that three of more sample 
  75. ;       populations have the same mean of distribution is rejected if two or
  76. ;       more populations differ with statistical significance. 
  77. ;
  78. ; REFERENCE:
  79. ;       PROBABILITY and STATISTICS for ENGINEERS and SCIENTISTS (3rd edition)
  80. ;       Ronald E. Walpole & Raymond H. Myers
  81. ;       ISBN 0-02-424170-9
  82. ;
  83. ; MODIFICATION HISTORY:
  84. ;       Written by:  GGS, RSI, September 1994
  85. ;-
  86.  
  87. pro crank, w, s
  88.   ;Replace elements of the sorted array "w" by their rank.
  89.   ;Identical observations ("ties") are ranked according to their means.
  90.   ;s = f^3 - f (f is the number of elements in identical observations.)
  91.   n = n_elements(w)
  92.   w = [0.0, w]  ;operate on elements w[1], ... , w[n] of the shifted
  93.                 ;n+1 element float array [w].
  94.   s = 0.0
  95.   j = 1
  96.   while(j lt n) do begin
  97.     if(w[j+1] ne w[j]) then begin
  98.       w[j] = j
  99.       j = j+1
  100.     endif else begin
  101.       for jt = j+1, n do $
  102.         if (w[jt] ne w[j]) then goto, case2
  103.       jt = n + 1
  104.       case2:
  105.       rank = 0.5 * (j + jt - 1)
  106.       for ji = j, jt-1 do $
  107.         w[ji] = rank
  108.       t = jt - j
  109.       s = s + t^3 - t
  110.       j = jt
  111.     endelse
  112.   endwhile
  113.   if(j eq n) then w[n] = n
  114.   w = w[1:*]
  115. end
  116.  
  117. function kw_test, x, df = df, missing = missing
  118.  
  119.   on_error, 2
  120.  
  121.   sx = size(x)
  122.   if sx[0] + 1 lt 3 then $
  123.     message, 'x must be a two-dimensional array with 3 or more columns.'
  124.  
  125.   ;Redimension x as a column vector.
  126.   xx = reform(transpose(x), 1, sx[4])
  127.  
  128.   if keyword_set(missing) eq 0 then begin
  129.     ;Equal length samples.
  130.     ixx = sort(xx) ;Sort and rank the combined vector.
  131.     xx = xx[ixx]
  132.     crank, xx
  133.     xx = xx[sort(ixx)]
  134.     crs = total(reform(xx, sx[2], sx[1]), 1) ;Compute the individual column 
  135.                                              ;rank sums by reforming xx.
  136.     n = sx[4]
  137.     h = 12.0/(n*(n+1))*total(crs^2 / sx[2]) - 3*(n+1)
  138.     df = sx[1] - 1
  139.     prob = 1 - chisqr_pdf(h, df)
  140.     return, [h, prob]
  141.   endif else begin
  142.     ;Unequal length samples.
  143.     ixx = where(xx ne missing, nxx) ;Eliminate missing data.
  144.     if nxx ne sx[4] then xx = xx[ixx]
  145.     ixx = sort(xx) ;Sort and rank the combined vector.
  146.     xx = xx[ixx]
  147.     crank, xx
  148.     xx = xx[sort(ixx)]
  149.     mvi = lonarr(sx[1])
  150.     if sx[3] eq 5 then crs = dblarr(sx[1]) else $
  151.       crs = fltarr(sx[1])
  152.     top = 0L
  153.     for k = 0, sx[1]-1 do begin 
  154.       imiss = where(x[k,*] eq missing, nmiss)
  155.       mvi[k] = nmiss ;Number of missing values per column vector.
  156.       bot = top + sx[2] - mvi[k] - 1L     
  157.       crs[k] = total(xx[top:bot]) ;Column rank sum.
  158.       top = bot + 1L
  159.     endfor
  160.     h = 12.0/(nxx*(nxx+1))*total(crs^2 / (sx[2]-mvi)) - 3*(nxx+1)
  161.     df = sx[1] - 1
  162.     prob = 1 - chisqr_pdf(h, df)
  163.     return, [h, prob]
  164.   endelse
  165.  
  166. end
  167.